1 Bug Reports

Please file issues here for any bugs or technical difficulties.

2 R Installation

R (version 3.4.3 or higher) is required: see https://www.r-project.org

Install packages:

source("http://bioconductor.org/biocLite.R")
biocLite(c("supraHex","devtools","dnet","XGR","dplyr","GenomicRanges"))
devtools::install_github(c("hfang-bristol/supraHex","hfang-bristol/XGR"))

IMPORTANT: Please check the package ‘GenomicRanges’ (version 1.30.0 or higher is required). If not, please obtain/update it from Bioconductor (https://www.bioconductor.org/install/)

3 Cubism (2D regular maps)

A collection of spatially ordered maps supported for the self-organised representation and comparison: (from top-left to bottom-right) a supra-hexagonal map, a ring map, a diamond map, a trefoil map, a butterfly map, an hourglass map, a ladder map, a bridge map, and a triangle map.

library(supraHex)
library(ggplot2)
colormap <- 'lightyellow-lightblue'
color <- 'black'

# For "supraHex" shape itself
sHex <- sHexGridVariant(r=6, shape="suprahex")
df_polygon <- sHexPolygon(sHex)
df_coord <- data.frame(sHex$coord, index=1:nrow(sHex$coord))
gp_suprahex <- ggplot(data=df_polygon, aes(x,y,group=index)) + geom_polygon(aes(fill=factor(stepCentroid%%2))) + coord_fixed(ratio=1) + theme_void() + theme(legend.position="none") + scale_fill_manual(values=visColormap(colormap)(2)) + geom_text(data=df_coord, aes(x,y,label=index), color=color, size=2.5) + theme(plot.margin=unit(rep(0,4),rep("lines",4)))

# For "triangle" shape
sHex <- sHexGridVariant(r=6, shape="triangle")
df_polygon <- sHexPolygon(sHex)
df_coord <- data.frame(sHex$coord, index=1:nrow(sHex$coord))
gp_triangle <- ggplot(data=df_polygon, aes(x,y,group=index)) + geom_polygon(aes(fill=factor(stepCentroid%%2))) + coord_fixed(ratio=1) + theme_void() + theme(legend.position="none") + scale_fill_manual(values=visColormap(colormap)(2)) + geom_text(data=df_coord, aes(x,y,label=index), color=color, size=2.5) + theme(plot.margin=unit(rep(0,4),rep("lines",4)))

# For "diamond" shape
sHex <- sHexGridVariant(r=6, shape="diamond")
df_polygon <- sHexPolygon(sHex)
df_coord <- data.frame(sHex$coord, index=1:nrow(sHex$coord))
gp_diamond <- ggplot(data=df_polygon, aes(x,y,group=index)) + geom_polygon(aes(fill=factor(stepCentroid%%2))) + coord_fixed(ratio=1) + theme_void() + theme(legend.position="none") + scale_fill_manual(values=visColormap(colormap)(2)) + geom_text(data=df_coord, aes(x,y,label=index), color=color, size=2.5) + theme(plot.margin=unit(rep(0,4),rep("lines",4)))

# For "hourglass" shape
sHex <- sHexGridVariant(r=6, shape="hourglass")
df_polygon <- sHexPolygon(sHex)
df_coord <- data.frame(sHex$coord, index=1:nrow(sHex$coord))
gp_hourglass <- ggplot(data=df_polygon, aes(x,y,group=index)) + geom_polygon(aes(fill=factor(stepCentroid%%2))) + coord_fixed(ratio=1) + theme_void() + theme(legend.position="none") + scale_fill_manual(values=visColormap(colormap)(2)) + geom_text(data=df_coord, aes(x,y,label=index), color=color, size=2.5) + theme(plot.margin=unit(rep(0,4),rep("lines",4)))

# For "trefoil" shape
sHex <- sHexGridVariant(r=6, shape="trefoil")
df_polygon <- sHexPolygon(sHex)
df_coord <- data.frame(sHex$coord, index=1:nrow(sHex$coord))
gp_trefoil <- ggplot(data=df_polygon, aes(x,y,group=index)) + geom_polygon(aes(fill=factor(stepCentroid%%2))) + coord_fixed(ratio=1) + theme_void() + theme(legend.position="none") + scale_fill_manual(values=visColormap(colormap)(2)) + geom_text(data=df_coord, aes(x,y,label=index), color=color, size=2.5) + theme(plot.margin=unit(rep(0,4),rep("lines",4)))

# For "ladder" shape
sHex <- sHexGridVariant(r=6, shape="ladder")
df_polygon <- sHexPolygon(sHex)
df_coord <- data.frame(sHex$coord, index=1:nrow(sHex$coord))
gp_ladder <- ggplot(data=df_polygon, aes(x,y,group=index)) + geom_polygon(aes(fill=factor(stepCentroid%%2))) + coord_fixed(ratio=1) + theme_void() + theme(legend.position="none") + scale_fill_manual(values=visColormap(colormap)(2)) + geom_text(data=df_coord, aes(x,y,label=index), color=color, size=2.5) + theme(plot.margin=unit(rep(0,4),rep("lines",4)))

# For "butterfly" shape
sHex <- sHexGridVariant(r=6, shape="butterfly")
df_polygon <- sHexPolygon(sHex)
df_coord <- data.frame(sHex$coord, index=1:nrow(sHex$coord))
gp_butterfly <- ggplot(data=df_polygon, aes(x,y,group=index)) + geom_polygon(aes(fill=factor(stepCentroid%%2))) + coord_fixed(ratio=1) + theme_void() + theme(legend.position="none") + scale_fill_manual(values=visColormap(colormap)(2)) + geom_text(data=df_coord, aes(x,y,label=index), color=color, size=2.5) + theme(plot.margin=unit(rep(0,4),rep("lines",4)))

# For "ring" shape
sHex <- sHexGridVariant(r=6, shape="ring")
df_polygon <- sHexPolygon(sHex)
df_coord <- data.frame(sHex$coord, index=1:nrow(sHex$coord))
gp_ring <- ggplot(data=df_polygon, aes(x,y,group=index)) + geom_polygon(aes(fill=factor(stepCentroid%%2))) + coord_fixed(ratio=1) + theme_void() + theme(legend.position="none") + scale_fill_manual(values=visColormap(colormap)(2)) + geom_text(data=df_coord, aes(x,y,label=index), color=color, size=2.5) + theme(plot.margin=unit(rep(0,4),rep("lines",4)))

# For "bridge" shape
sHex <- sHexGridVariant(r=6, shape="bridge")
df_polygon <- sHexPolygon(sHex)
df_coord <- data.frame(sHex$coord, index=1:nrow(sHex$coord))
gp_bridge <- ggplot(data=df_polygon, aes(x,y,group=index)) + geom_polygon(aes(fill=factor(stepCentroid%%2))) + coord_fixed(ratio=1) + theme_void() + theme(legend.position="none") + scale_fill_manual(values=visColormap(colormap)(2)) + geom_text(data=df_coord, aes(x,y,label=index), color=color, size=2.5) + theme(plot.margin=unit(rep(0,4),rep("lines",4)))

gridExtra::grid.arrange(grobs=list(gp_suprahex,gp_ring,gp_diamond,gp_trefoil,gp_butterfly,gp_hourglass,gp_ladder,gp_bridge,gp_triangle), layout_matrix=rbind(c(1,1,2,2,3),c(1,1,2,2,3),c(4,4,5,5,6),c(4,4,5,5,6),c(7,7,8,8,9)), nrow=5, ncol=5)

4 Showcase-I

In this showcase, we analyse TF and DHS datasets in K562, and compare TF datasets between K562 and Gm12878.

First of all, load the packages:

library(supraHex)
library(XGR)

# Specify the local location of built-in RData
RData.location <- "http://galahad.well.ox.ac.uk/bigdata"

Then, import TF and DHS datasets:

# extract TFBSs, stored in a list of GR objects
ENCODE_TFBS_ClusteredV3_CellTypes <- xRDataLoader(RData.customised='ENCODE_TFBS_ClusteredV3_CellTypes', RData.location=RData.location)
## common to Gm12878 and K562
anno_gr_Gm12878 <- ENCODE_TFBS_ClusteredV3_CellTypes$GM12878
anno_gr_K562 <- ENCODE_TFBS_ClusteredV3_CellTypes$K562
names_common <- sort(intersect(names(anno_gr_Gm12878), names(anno_gr_K562)))
ind <- match(names_common, names(anno_gr_Gm12878))
TF_gr_Gm12878 <- anno_gr_Gm12878[ind]
ind <- match(names_common, names(anno_gr_K562))
TF_gr_K562 <- anno_gr_K562[ind]

# extract DHS, stored in a list of GR objects
ENCODE_DNaseI_ClusteredV3_CellTypes <- xRDataLoader(RData.customised='ENCODE_DNaseI_ClusteredV3_CellTypes', RData.location=RData.location)
## in K562
ind <- grep("K562", names(ENCODE_DNaseI_ClusteredV3_CellTypes))
DHS_gr_K562 <- ENCODE_DNaseI_ClusteredV3_CellTypes[ind]

4.1 Learned TF map in K562

Calculate gene-centric TF association scores:

anno_gr <- TF_gr_K562
res_ls <- lapply(1:length(anno_gr), function(i){
    message(sprintf("Analysing %d '%s' (%s)", i, names(anno_gr)[i], as.character(Sys.time())), appendLF=T)
    df_xGenes <- xGR2xGenes(anno_gr[[i]], format="GRanges", crosslink="genehancer", cdf.function="empirical", scoring=T, scoring.scheme="max", RData.location=RData.location, verbose=F)
    data.frame(df_xGenes, TF=rep(names(anno_gr)[i],nrow(df_xGenes)), stringsAsFactors=F)
})
df <- do.call(rbind, res_ls)
G2TF_K562 <- as.matrix(xSparseMatrix(df[,c("Gene","TF","Score")]))

Justify the use of the ladder-shaped map:

data <- G2TF_K562
df <- as.data.frame(stats::prcomp(data)$x[,1:2])
gp <- ggplot(df, aes(x=PC1, y=PC2)) 
gp <- gp + theme_bw() + theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank())
## 2d shape plot
myColor <- xColormap(colormap="spectral")(11)
gp <- gp + stat_binhex(bins=48) + scale_fill_gradientn(colours=myColor,trans="sqrt") + geom_density2d(colour="grey")
print(gp)

Learn the ladder-shaped map using cross-TF gene scores:

data <- G2TF_K562
sMap <- sPipeline(data, shape="ladder", finetuneSustain=T)
#sWriteData(sMap=sMap, data=data, filename="output.G2TF_K562.txt", keep.data=T)

Build neighbour-joining TF tree:

D <- t(sMap$codebook)
set.seed(825)
tree_bs <- visTreeBootstrap(D, num.bootstrap=1000, nodelabels.arg=list(frame="circle",cex=0.35,col="black"), plot.phylo.arg=list(type=c("phylogram","cladogram","fan","radial")[1],cex=0.6,node.pos=1))

Visualise TF map (reordered by tree):

tree_bs_ind <- match(tree_bs$tip.label, rownames(D))
xMap <- sMap
xMap$codebook <- xMap$codebook[,tree_bs_ind]
colormap <- xColormap("RdYlBu")
visHexMulComp(xMap, which.components=51:1, rect.grid=c(7,8), title.rotate=0, title.xy=c(0.3,1.05), colormap=colormap, zlim=c(0,1), gp=grid::gpar(cex=0.6), newpage=F)

4.2 Fused DHS map in K562

Calculate gene-centric DHS association scores:

anno_gr <- DHS_gr_K562
res_ls <- lapply(1:length(anno_gr), function(i){
    message(sprintf("Analysing %d '%s' (%s)", i, names(anno_gr)[i], as.character(Sys.time())), appendLF=T)
    df_xGenes <- xGR2xGenes(anno_gr[[i]], format="GRanges", crosslink="genehancer", cdf.function="empirical", scoring=T, scoring.scheme="max", RData.location=RData.location, verbose=T)
    data.frame(df_xGenes, Cell=rep(names(anno_gr)[i],nrow(df_xGenes)), stringsAsFactors=F)
})
df <- do.call(rbind, res_ls)
G2DHS <- as.matrix(xSparseMatrix(df[,c("Gene","Cell","Score")]))
#output <- data.frame(Gene=rownames(G2DHS), G2DHS, stringsAsFactors=F)
#write.table(output,file='output.G2DHS.txt',sep='\t',row.names=F,quote=F)

The DHS map is obtained by overlaying/fusing the DHS data onto the learned TF map in K562:

ind <- match(rownames(G2TF_K562), rownames(G2DHS))
additional <- G2DHS[ind,]
additional[is.na(additional)] <- 0
additional <- matrix(additional, ncol=1)
rownames(additional) <- rownames(G2TF_K562)
colnames(additional) <- 'DHS'
sOverlay_DHS <- sMapOverlay(sMap, data=G2TF_K562, additional=additional)

Visualise DHS map:

colormap <- xColormap("RdYlBu")
visHexMulComp(sOverlay_DHS, rect.grid=c(1,1), title.rotate=0, colormap=colormap, zlim=c(0,1), gp=grid::gpar(cex=1), newpage=F)

Rankings of 51 TFs in terms of similarity (Pearson correlation) to the DHS map in K562:

M_K562 <- sMap$codebook
M_DHS <- sOverlay_DHS$codebook
y <- M_DHS[,1]
names(y) <- 1:length(y)
ls_df <- lapply(1:ncol(M_K562), function(i){
    x <- data.frame(name=1:nrow(M_K562), M_K562[,i], stringsAsFactors=F)
    ls_res <- xCorrelation(x, y, method="pearson", p.type="nominal")
    df_summary <- ls_res$df_summary
    data.frame(TF=colnames(M_K562)[i], cor=df_summary$cor, stringsAsFactors=F)
})
df <- do.call(rbind, ls_df)
df$rank <- rank(-1*df$cor,ties.method="first")
df <- df %>% dplyr::arrange(rank)

## polar plot
gp_DHS <- xPolarDot(df, parallel=T, signature=F)
gp_DHS + labs(y="Pearson correlation") + ylim(0,1)

4.3 Fused TF map in Gm12878

Calculate gene-centric TF association scores:

anno_gr <- TF_gr_Gm12878
res_ls <- lapply(1:length(anno_gr), function(i){
    message(sprintf("Analysing %d '%s' (%s)", i, names(anno_gr)[i], as.character(Sys.time())), appendLF=T)
    df_xGenes <- xGR2xGenes(anno_gr[[i]], format="GRanges", crosslink="genehancer", cdf.function="empirical", scoring=T, scoring.scheme="max", RData.location=RData.location, verbose=F)
    data.frame(df_xGenes, TF=rep(names(anno_gr)[i],nrow(df_xGenes)), stringsAsFactors=F)
})
df <- do.call(rbind, res_ls)
G2TF_Gm12878 <- as.matrix(xSparseMatrix(df[,c("Gene","TF","Score")]))
#output <- data.frame(Gene=rownames(G2TF_Gm12878), G2TF_Gm12878, stringsAsFactors=F)
#write.table(output,file='output.G2TF_Gm12878.txt',sep='\t',row.names=F,quote=F)

The TF map is obtained by overlaying/fusing the TF data in Gm12878 onto the learned TF map in K562:

ind <- match(rownames(G2TF_K562), rownames(G2TF_Gm12878))
additional <- G2TF_Gm12878[ind,]
additional[is.na(additional)] <- 0
sOverlay_Gm12878 <- sMapOverlay(sMap, data=G2TF_K562, additional=additional)

Rankings of 51 TFs in terms of similarity (Pearson correlation) in two cell lines (K562 and Gm12878):

M_K562 <- sMap$codebook
M_Gm12878 <- sOverlay_Gm12878$codebook
ls_df <- lapply(1:ncol(M_K562), function(i){
    x <- data.frame(name=1:nrow(M_K562), M_K562[,i], stringsAsFactors=F)
    y <- M_Gm12878[,i]
    names(y) <- 1:nrow(M_Gm12878)
    ls_res <- xCorrelation(x, y, method="pearson", p.type="nominal")
    df_summary <- ls_res$df_summary
    data.frame(TF=colnames(M_K562)[i], cor=df_summary$cor, stringsAsFactors=F)
})
df <- do.call(rbind, ls_df)
df$rank <- rank(-1*df$cor,ties.method="first")
df <- df %>% dplyr::arrange(rank)

## polar plot
gp_twin <- xPolarDot(df, parallel=T, signature=F)
gp_twin + labs(y="Pearson correlation") + ylim(0,1)

Visualise TF map (reordered by correlation):

## yMap
yMap <- sMap
ind <- match(df$TF, colnames(yMap$codebook))
yMap$codebook <- yMap$codebook[,ind]
## zMap
zMap <- sOverlay_Gm12878
ind <- match(df$TF, colnames(zMap$codebook))
zMap$codebook <- zMap$codebook[,ind]
## xMap
xMap <- sMap
nr <- nrow(xMap$codebook)
nc <- ncol(xMap$codebook)
matrix.combined <- matrix(NA, nrow=nr, ncol=nc*2)
matrix.combined[, seq(1,nc*2,2)] <- yMap$codebook
matrix.combined[, seq(2,nc*2,2)] <- zMap$codebook
xMap$codebook <- matrix.combined
colnames(xMap$codebook) <- rep(colnames(yMap$codebook), each=2)
## visualisation
colormap <- xColormap("RdYlBu")
visHexMulComp(xMap, which.components=1:(nc*2), rect.grid=c(13,8), title.rotate=0, title.xy=c(0.3,1), colormap=colormap, zlim=c(0,1), gp=grid::gpar(cex=0.4), newpage=F)

Perform enrichment analysis for TFs:

## TFs (high correlation)
data_high <- subset(df, rank<=25)$TF
## TFs (low correlation)
data_low <- subset(df, rank>=27)$TF

## enrichment analysis
data <- list(high=data_high, low=data_low)
ls_eTerm <- xEnricherGenesAdv(data, ontologies="SF", ontology.algorithm="none", RData.location=RData.location)

## visualisation
gp <- xEnrichForest(ls_eTerm, FDR.cutoff=0.01, zlim=c(2,10), signature=F)
gp + theme(axis.text.y=element_text(size=10))

Enriched protein domains for TFs

ls_eTerm$df

5 Showcase-II

In this showcase, we analyse TF and CRISPR datasets in K562.

First of all, load the packages:

library(supraHex)
library(XGR)

# optionally, specify the local location of built-in RData
# by default:
RData.location <- "http://galahad.well.ox.ac.uk/bigdata_dev"

Then, import TF and CRISPR datasets in K562:

# extract TF, stored in a list of GR objects
ENCODE_TFBS_ClusteredV3_CellTypes <- xRDataLoader(RData.customised='ENCODE_TFBS_ClusteredV3_CellTypes', RData.location=RData.location)
# calculate gene-centric TF association scores
anno_gr <- ENCODE_TFBS_ClusteredV3_CellTypes$K562
res_ls <- lapply(1:length(anno_gr), function(i){
    message(sprintf("Analysing %d '%s' (%s)", i, names(anno_gr)[i], as.character(Sys.time())), appendLF=T)
    df_xGenes <- xGR2xGenes(anno_gr[[i]], format="GRanges", crosslink="genehancer", cdf.function="empirical", scoring=T, scoring.scheme="max", RData.location=RData.location, verbose=F)
    data.frame(df_xGenes, TF=rep(names(anno_gr)[i],nrow(df_xGenes)), stringsAsFactors=F)
})
df <- do.call(rbind, res_ls)
mat_crispr <- as.matrix(xSparseMatrix(df[,c("Gene","TF","Score")]))

# import CRISPR
data.file <- file.path(RData.location,'TW_Science_CRISPR.txt')
data <- read.delim(data.file, header=TRUE, stringsAsFactors=FALSE)
data <- subset(data, K562.adjusted.p.value<0.05)

# extract genes commons to CRISPR and TF datasets in K562
ind <- match(data$Gene, rownames(mat_crispr))
G2TF_K562_crispr <- mat_crispr[ind[!is.na(ind)],]
G2CRISPR_K562 <- data.frame(CRISPR=data$K562.CS[!is.na(ind)], stringsAsFactors=F)
rownames(G2CRISPR_K562) <- data$Gene[!is.na(ind)]

5.1 Learned TF map in K562

Justify the use of the triangle-shaped map:

data <- G2TF_K562_crispr
df <- as.data.frame(stats::prcomp(data)$x[,1:2])
gp <- ggplot(df, aes(x=PC1, y=PC2)) 
gp <- gp + theme_bw() + theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank())
## 2d shape plot
myColor <- xColormap(colormap="spectral")(11)
gp <- gp + stat_binhex(bins=24) + scale_fill_gradientn(colours=myColor,trans="sqrt") + geom_density2d(colour="grey")
print(gp)

Fuse the triangle-shaped map using cross-TF gene scores:

data <- G2TF_K562_crispr
sM <- sPipeline(data, shape="triangle", finetuneSustain=T, scale=3)

5.2 Fused CRISPR map in K562

The CRISPR map is obtained by overlaying/fusing the CRISPR data onto the learned TF map in K562:

sOverlay_CRISPR <- sMapOverlay(sM, data=G2TF_K562_crispr, additional=G2CRISPR_K562)

Visualise CRISPR map:

colormap <- xColormap("ggplot2.bottom")
visHexMulComp(sOverlay_CRISPR, rect.grid=c(1,1), title.rotate=0, title.xy=c(0.4,1), colormap=colormap, zlim=c(-2.5,-1.5), gp=grid::gpar(cex=1), newpage=F)

Rankings of 100 TFs in terms of similarity (Pearson correlation) to the CRISPR map in K562:

M_K562 <- sM$codebook
M_CRISPR <- sOverlay_CRISPR$codebook
y <- M_CRISPR[,1]
names(y) <- 1:length(y)
ls_df <- lapply(1:ncol(M_K562), function(i){
    x <- data.frame(name=1:nrow(M_K562), M_K562[,i], stringsAsFactors=F)
    ls_res <- xCorrelation(x, y, method="pearson", p.type="nominal")
    df_summary <- ls_res$df_summary
    data.frame(TF=colnames(M_K562)[i], cor=df_summary$cor, stringsAsFactors=F)
})
df <- do.call(rbind, ls_df)
df$rank <- rank(-1*df$cor,ties.method="first")
df <- df %>% dplyr::arrange(rank)

## polar plot
gp_CRISPR <- xPolarDot(df, size=1.5, parallel=T, signature=F)
gp_CRISPR + labs(y="Pearson correlation") + scale_y_reverse(limits=c(0,-1))

Visualise TF map (reordered by correlation):

xMap <- sM
ind <- match(df$TF, colnames(xMap$codebook))
xMap$codebook <- xMap$codebook[,ind]
colormap <- xColormap("RdYlBu")
visHexMulComp(xMap, which.components=100:1, rect.grid=c(10,10), title.rotate=0, title.xy=c(0.3,1.05), colormap=colormap, zlim=c(0,1), gp=grid::gpar(cex=0.6), newpage=F)

5.3 Gene clusters

Obtain gene clusters based on TF map:

xM <- sM
ind <- match(df$TF, colnames(xM$codebook))
xM$codebook <- xM$codebook[,ind]
#seed <- sDmatMinima(xM)
sBase <- sDmatCluster(xM, which_neigh=1)
colormap <- xColormap("spectral")
visDmatCluster(xM, sBase, colormap=colormap, newpage=F)

Calculate CRISPR scores per gene cluster:

colormap <- xColormap("grey-yellow-red")
output <- visDmatHeatmap(xM, data=G2TF_K562_crispr[,ind], sBase, base.color=xColormap("spectral"), base.separated.arg=list(col="white",lty=1,lwd=1), base.legend.location="bottomleft", colormap=colormap, zlim=c(0,1), KeyValueName="TF-gene association", labRow=NA, srtCol=90, cexCol=0.6, keep.data=T)

## Polar barplot
ind <- match(output$ID, rownames(G2CRISPR_K562))
output$val <- G2CRISPR_K562[ind,]
#write.table(output,file='output.CRISPR.txt',sep='\t',row.names=F)
ls_tmp <- split(x=output$val, f=output$Cluster_base)
ls_df <- lapply(1:length(ls_tmp), function(i){
    x <- ls_tmp[[i]]
    y <- median(x)
    data.frame(Cluster=paste0('C',names(ls_tmp)[i]), Val=y, Num=length(x), stringsAsFactors=F)
})
df <- do.call(rbind, ls_df)
gp <- xPolarBar(df, colormap='spectral', size.name=10, size.value=2.5, parallel=F, signature=F)
gp

Perform enrichment analysis for clustered genes:

ls_cluster <- split(x=output$ID, f=output$Cluster_base)
names(ls_cluster) <- paste0('C', names(ls_cluster))
background <- output$ID
ontologies <- "REACTOME"
ls_res_reactome <- xEnricherGenesAdv(list_vec=ls_cluster, background=background, ontologies=ontologies, size.range=c(20,500), test="fisher", min.overlap=10, p.tail="one-tail", RData.location=RData.location, plot=T, fdr.cutoff=0.05, displayBy="zscore")
#ls_res_reactome$gp + theme(legend.title=element_text(size=8),axis.text.x=element_text(size=4))
gp <- xEnrichForest(ls_res_reactome, CI.one=F, FDR.cutoff=0.05, signature=F, drop=T, colormap="yellow-red")
gp